[BZOJ1013]-球形空间产生器

传送门

题意略。

对于每个点，他们到球心的距离都是相等的。

根据欧几里得距离，此时只需求出一点 $(x_1, x_2, … x_n)$，使得：

$$\sum\limits_{j = 1}^n (a_{i, j} - x_j) ^ 2 = C$$

其中 $C$ 为常数。该方程组是由 $n + 1$ 个 $n$ 元二次方程组成的，不是线性方程组（线性方程组要求均为一次），但我们可以通过相邻两个方程作差，把它变成 $n$ 个 $n$ 元一次方程，同时消去常数 $C$：

$$\sum\limits_{j = 1}^n (a_{i, j}^2 - a_{i + 1, j}^2 - 2x_j(a_{i, j} - a_{i + 1, j})) = 0\ \ \ \ (i = 1, 2, ... , n)$$

把变量放在左边，常数放在右边：

$$\sum\limits_{j = 1}^n 2(a_{i, j} - a_{i + 1, j})x_j = \sum\limits_{j = 1}^n (a_{i, j}^2 - a_{i + 1, j}^2)\ \ \ \ (i = 1, 2, ... , n)$$

这就是一个线性方程组了，题目保证方程组有唯一解，我们直接对下面的增广矩阵进行高斯消元，变为简化阶梯形矩阵。

$$\begin{bmatrix} 2(a_{1, 1} - a_{2, 1}) & 2(a_{1, 2} - a_{2, 2}) & \cdots & 2(a_{1, n} - a_{2, n}) & \sum_{j = 1}^n(a_{1, j}^2 - a_{2, j}^2) \\ 2(a_{2, 1} - a_{3, 1}) & 2(a_{2, 2} - a_{3, 2}) & \cdots & 2(a_{2, n} - a_{3, n}) & \sum_{j = 1}^n(a_{2, j}^2 - a_{3, j}^2) \\ \vdots & \vdots & \ddots & \vdots & \vdots \\ 2(a_{n, 1} - a_{n + 1, 1}) & 2(a_{n, 2} - a_{n + 1, 2}) & \cdots & 2(a_{n, n} - a_{n + 1, n}) & \sum_{j = 1}^n(a_{n, j}^2 - a_{n + 1, j}^2) \\ \end{bmatrix}$$

code:

#include <bits/stdc++.h>
using namespace std;

double a[20][20], b[20], c[20][20];
int n;

int main() {
    scanf("%d", &n);
    for (int i = 1; i <= n + 1; i++)
        for (int j = 1; j <= n; j++) scanf("%lf", &a[i][j]);
    for (int i = 1; i <= n; i++)
        for (int j = 1; j <= n; j++) {
            c[i][j] = 2 * (a[i][j] - a[i + 1][j]);
            b[i] += a[i][j] * a[i][j] - a[i + 1][j] * a[i + 1][j];
        }
    for (int i = 1; i <= n; i++) {
        for (int j = i; j <= n; j++) {
            if (fabs(c[j][i]) > 1e-8) {
                for (int k = 1; k <= n; k++) swap(c[i][k], c[j][k]);
                swap(b[i], b[j]);
                break;
            }
        }
        for (int j = 1; j <= n; j++) {
            if (i == j) continue;
            double rate = c[j][i] / c[i][i];
            for (int k = i; k <= n; k++) c[j][k] -= c[i][k] * rate;
            b[j] -= b[i] * rate;
        }
    }
    for (int i = 1; i <= n; i++) printf("%.3lf ", b[i] / c[i][i]);
    cout << endl;
    return 0;
}

高斯消元不怎么好写，重点是怎么把题意转化为式子，还有推式子的时候要勤快！